####################################################
####          Replication Code for              ####
####     Bayesian versus maximum likelihood     ####
####      estimation of treatment effects       ####
####          in bivariate probit               ####  
####       instrumentalvariable models          ####    
####################################################
####                                            ####  
####                3/22/2017                   ####         
####     this file runs rcode to create Figures ####  
####    and creates a log file                  ####  
####################################################

#install.packages(c("abind", "ggplot2", "gridExtra", "grDevices", "grid", "ggsn"))
### create log file
sink("log_SimulationResults.txt")



library(abind)
library(ggplot2)
library(gridExtra)
library(grDevices)
library(grid)
library(rstan)
library(xtable)
library(plyr)
library(pbivnorm)
library(haven)
library(ggsn)
sessionInfo()
### before running set working directory to folder with files


### functions for LATE, ATE, MAE. MSE, 
source("HelperFunctions.R", echo = T)

## recreate sim data
N_seq = c(50, 100, 250, 500)
cor = c(0.25,0.5,0.75)
IV_strength = c(0.5, 1, 1.5, 2)
treatment = c(0, 1, 2, 3, 4)
simulation_para = data.frame(do.call(rbind,replicate(500,as.matrix(expand.grid(N_seq, cor, IV_strength, treatment)),simplify = F)))
names(simulation_para) = c("N","cor","IV_strength","treatment")
simulation_para = simulation_para[order(simulation_para$N, simulation_para$cor, simulation_para$IV_strength, simulation_para$treatment),]


### create plots
source("Figures_MainText_AppendixG.R", echo = T) # good
source("Figures_AppendixE.R", echo = T) 
source("Figures_AppendixF.R", echo = T)

sink()
